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low temperatures, a regime notoriously difficulty to simulate because the large free- 
energy barriers. The good results found (when compared with other well established 
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phase transitions, a possibility still not well explored in the literature. 
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I. INTRODUCTION 



Many statistical systems are difficult to "probe" since their phase spaces display com- 
plicated landscapes full of energetic valleys and hills^. In such case, a non-representative 
sampling of the microstates, e.g., due to uneven visits to the different domains^, can lead to 
metastability and broken ergodicity, with a consequent non-convergence to equilibrium and 
poor estimates for the thermodynamic quantities 4,5 . This is exactly the situation found in 
first-order phase transitions 6 , where the free-energy minima are separated by large barriers, 
and simple one-flip Metropolis approaches are unable to solve the problem. 

Thus, different methods - aimed to guarantee ergodic simulations in the context described 
above - have been proposed^ - —. Among the so called enhanced sampling algorithms, a 
particularly important one due to its simplicity and generality is the simulated tempering 
(ST)^£ (closely related to the also relevant parallel (or replica) tempering, PT, approach^). 
Here, tempering means that along the simulations the system can undergo temperature 
changes. Consider we shall analyze a system at T = T\. Besides an usual Monte Carlo (MC) 
prescription, in the ST algorithm T is also treated as a dynamical variable: the temperature 
can assume distinct values from a set {T% < T 2 < . . . < T/v}, switching from time to time 
according to established accepting probabilities p's. Of course, during the simulations the 
relevant averages necessary to obtain the sought thermodynamic quantities are performed 
only when the system is at T\. The central idea is that temporary evolutions at higher T's 
strongly facilitates the crossing of the free-energy barriers, then allowing uniform visits to 
the multiple regions of a fragmented phase-space^. 

A fundamental ingredient in the ST is the definition of the accepting probabilities^; 
functions of the energy, the different T's and appropriate weight factors g (see next Section 
as well as interesting discussions in Ref.—). By their turn, the exact expressions for the 
g's depend on the problem partition function Z, a quantity often difficult to calculate^ - —, 
demanding computationally time consuming methods^. Obviously, approximations for the 
g's can be used (e.g., as in Ref.—). But so the ST efficiency can be dramatically hindered, 
once the weights themselves can create certain bias for the sampling. For instance, they may 
lead to unbalanced sorted temperatures, given rise to a non-ergodic visit of the phase-space if 
the probabilities to pick the lower T n 's are too high. On the other hand, if the more frequent 
T's are the greater ones, the generated set of microstates cannot properly characterize the 

2 



system features at T = Ti, the actual temperature of interest. Such technical aspects 
associated to the g's are even more delicate for first order phase transitions, a situation 
rarely analyzed with the ST->22. 

Given the previous comments, our purpose in this contribution is twofold. First, to 
present a numerically simple - yet quite accurate - procedure to obtain the weight factors 
g, thus making the ST algorithm easier to implement, e.g., by avoiding involving recursive 
protocols^. This is accomplished with a method proposed in Ref.-^, where Z is calculated 
directly from the transfer matrix largest eigenvalue (A*- )), straightforward to compute from 
Monte Carlo simulations. The key point is that although gives the exact Z only at the 
thermodynamical limit, i.e., for infinite systems, the convergence is very fast. So, even for 
a relatively small system (as will be illustrated in the examples), any difference between its 
exact Z and that from can be neglected and for all practical reasons the approach leads 
to the correct g's. Second, to consider the ST for the already mentioned difficult case of 
first-order transitions, showing that the ST is also a helpful tool to address such regime, a 
possibility barely explored in the literature. 

The work is organized as the following. We review the ST algorithm and how to calculate 
g from the transfer matrix in Section II. Also in Section II we exemplify the procedure 
efficiency by computing the partition function for the Ising model, a system for which Z can 
be obtained exactly. In Section III and IV we compare the present with other well established 
methods to study first-order phase transitions, taking as case studies the Blume-Capel, 
Blume-Emery-Griffiths (BEG) and Bell-Lavis models. Finally, remarks and the conclusion 
are drawn in Section V. 

II. THE ST AND THE CALCULATION OF g 

As mentioned in the Introduction, the ST algorithm is generally implemented as two 
steps procedure, repeated a given number of times. First, at a temperature T n > (n' = 
1, 2, . . . , N), a standard Metropolis prescription is used to promote the transition a' — > a" 
with the probability P a '^ a » = min{l, exp[— f3 n >('H(a") — *H(<r'))]} (for o representing the 
system microscopic configurations). Second, an attempt to change the replica temperature 
(from T n i to T n ») is made according to 

Vn'^n" = min{l, exp[(/3 n / - p n »)K(a) + (g n >> ~ 9n'))}- (1) 
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In Eq. (PQ), H is the problem Hamiltonian, a the actual microscopic state, and n" = 1, . . . , N 
arbitrary, so non-adjacent temperatures changes are allowed. Moreover, = 1/(/cbT) with 
k B always set equal to 1 hereafter. 

We observe that p n '^n" is strongly dependent on the weights g's. Indeed, for an appro- 
priate sampling the evolution should uniformly visit all the established temperatures, which 
is the case when g n = (3 n f n for f n the free-energy per volume V at T n . Recalling the rela- 
tion P n f n = — \n[Z n ]/V, we see that the calculation of g is not a trivial task: neither the 
partition function nor the free-energy can be obtained directly from MC simulations since 
there are no thermodynamic quantities whose averages lead to Z and /. For this reason 
some alternative methods have been proposecU^^ 1 ^. Here we shall consider a rather simple 
numerical approach to compute frr^-, based on the transfer matrix T largest eigenvalue \(°\ 

Briefly, at the thermodynamic limit it holds true that (see details in the Appendix) 

Zn = (^) K , (2) 

with K ~ V l l d and d the spatial dimension. To obtain A^ is straightforward. In fact, 
suppose for defmiteness a 2D system with K layers of L sites each, so V = L x K. Next, 
consider the full Hamiltonian decomposed as 

k=K 

H= 'HjSk, Sk+i), (3) 
fc=i 

where = (cr ljfc , a 2; k, ■ ■ ■ , <?L,k) represents the state configuration of the fc-th layer. We 
further assume periodic boundary conditions, or Sk+i = Si. The transfer matrix is defined 
(Appendix) so that its elements reacr^ 

T(S k , S k+ x) = exp[-f3U(S k , S k+ x)]. (4) 

Then, as shown in the Appendix, we have (T(Sk, Sfc+i = Sfc) = T(Sk)) 

A (0) = (T(S k ))/(5 SkA+1 ). (5) 

This expression enables one to calculate the largest eigenvalue of T in terms of the 
averages (T(S k )) and (8 Sk ,s k+1 ), with S Sk ,s k+1 = 1 ($s k ,s k+1 = 0) if the layers S k and S k+ i 
are equal (different). We also should mention that (T(S k )) and {Ss k ,s k+1 ) can be evaluated 
from quite standard MC simulations. 
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As the final step, the weights follow from 



(6) 



A relevant issue is that even thought Eq. is exactly only for infinite size systems, if L 
(or K) is not too small the relation is extremely accurate. Hence, for any practical purpose 
Eq. gives the correct g, as we are going to illustrate in the next Sections. 

In this way, we can summarize the proposed approach as the following. First, with Eqs. 
(JSJl-flSD one evaluates the partition function and consequently the free-energy weights (at the 
temperatures {T n }) from usual MC simulations. Second, having the correct g's, one just 
implement the ST algorithm as previously explained. 

Finally, we comment on three important technical issues. The first is related to the query: 
how can standard MC simulations lead to good values for g around first-order phase tran- 
sition if then Metropolis algorithm usually yields unbalanced samplings? The answer relies 
on the fact that the free-energies are all equal at such regime of phases coexistence. Hence, 
even a biased microscopic sampling will result in accurate free-energies and consequently 
appropriate g's. On the other hand, the metastability typical of T's in the vicinity of the 
transition temperature can difficult the determination (with the necessary precision) of the 
function Z x T in such interval. Thus, derivatives of Z with respect to distinct parameters, 
representing different thermodynamic quantities, will give poor results. So, a second point is 
that the transfer matrix alone is not a reliable method to study first-order phase transitions. 
Lastly, we mention that in the limit of large systems and high temperatures (5s k ,s k+1 ) is 
small, since the probability for the configurations and S^+i to be the same is very low. 
As a consequence, one may get bad estimations for the weights. In this case, a possible way 
to circumvent the problem is to decrease the size L of each layer and to increase the number 
of layers, maintaining the volume V = L x K constant (see also the discussion in the Section 



A. An example: the partition function for the Ising model 

Just to verify how good is the Eq. (J2J) for finite size systems, we consider the Ising model, 
whose partition function can be calculated exactly. The Hamiltonian is 
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FIG. 1. For the Ising model, the partition function versus T calculated exactly** and from Eq. 
([2]) (symbols), with the parameters as in the text. The exact solutions, the lines, for L = 16 
(continuous) and L = 32 (dotted) are almost indistinguishable. 

where < i,j > denotes nearest-neighbors pairs i and j in a lattice of V = L d sites. At each 
site i, the spin variable assumes the values = ±1. J is the interaction energy and H is 
the magnetic field. For a square lattice (d = 2), T(Sk) yields 



Numerical calculations have been performed for H = (and in units of J). For L = 16 
and L = 32, Fig. [1] compares the exact partition function for a finite system (obtained 
from the solution in Ref.— ) with that calculated from Eq. (J2J). The agreement is indeed 
remarkable, indicating that even for relatively small systems, Z and therefore / are already 
very close to their values at the thermodynamic limit. 

Since the above system presents a very well known and simple first-order phase transition 
(for H = and T < T c ), we prefer to address such regime, our focus in this contribution, 
for other models in the following Sections. 

III. THE LATTICE-GAS MODEL WITH VACANCIES (BEG) 

The lattice-gas with vacancies (BEG) model is given by the Hamiltonian (er = 0, ±1) 



L 



T{S k ) = exp /3 ( ^2 J (1 + oi )k <Ti+i >k ) + H cri )k ) 



(8) 



i=i 
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for which the transfer matrix T(Sk) elements are 

T(S h ) = exp [fljr ((H + Jtn+u,) a hk + (J - D + K (1 + af +1}k )) a^)] . (10) 
1=1 

For the values of Kj J we are going to consider here and at low temperatures, if D is small, 
the system presents an ordered phase. When D increases, a gas phase takes place. These 
regimes are separated by a strong first-order phase transition at D = D*. For simplicity, 
hereafter all the quantities will be given in units of J. 



A. The Blume-Capel model 

For the Blume-Capel case, K = 0, accurate estimates for / (which directly gives the 
ST weights once g = (3 f) is available from the very efficient Wang-Landau method. We 
then compare in Fig. [2] the free-energy from our proposed procedure with that from Wang- 
Landau's^, considering K = 0, L = 32 and different .D's. As for the Ising model, again we 
see an excellent agreement (even for D = 1.965, the value corresponding to the tricritical 
point for T about 0.609). We also have explicit tested smaller L's (down to 12), always 
obtaining very good results. So, the exact g's for the ST are adequately (and easily) obtained 
from the transfer matrix approach. 

To study the Blume-Capel model around first-order phase transition, we perform finite 
size scaling analysis using the PT algorithm. We first recall that according to a rigorous 
theory of first-order phase transitions at low temperatures^, all thermodynamic quantities 
should scale with V. Then, for Q = Yn=i a h i n Fig. [3] we show the order parameter 
1 — (Q) /V an d the isothermal susceptibility xt = (^s^) _1 ((Q 2 ) — (Q) 2 )/V as functions 
of D for three system sizes L. The transition point can be estimated, for instance, from 
the peak position in the susceptibility or in the specific heat curves. Here we obtain D* 
through the location of the distinct L isotherms crossing*^-— . We find that for T = 0.50 
all the isotherms cross at D* = 1.9879(1), which closely agree with the values T = 0.499(3) 
and D* = 1.992(1) in Ref.— . In addition, the good collapse of the data in the insets, 
plotted as q x (D — D*) V and as Xt/V x (D — D*) V, illustrates the adequacy of the ST 
to locate transition points. Of course, other exact implementations of the ST, like that in 
Ref.—, would also solve the problem, however, by means of more complicated protocols, thus 
demanding longer computational times. 
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FIG. 2. For the Blume-Capel model, the free-energy versus T for different D obtained from the 
present approach (symbols) and from the Wang-Landau method (continuous lines). 
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FIG. 3. For the Blume-Capel model, the order parameter q and the isothermal susceptibility xt 
versus D for different system sizes at T = 0.5. The insets show the collapsed data, respectively, 
plotted as q x (D - D*) V and xt/V X (D - D*) V. 

B. The full, K ^ 0, BEG model 

Next, we briefly comment on the full BEG model, i.e., K ^ 0. Generally, approaches 
based on cluster algorithms are known to be very appropriate for lattice-gas systems^. Nev- 
ertheless, it is also a fact that first-order phase transitions for some special i^'s, like K = 3.3, 
can be a little trick to solve. So, modifications in the cluster method are necessary^. 
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FIG. 4. For the BEG model, with K/J = 3.3 and parameters as in the text, the probability 
distribution histogram of the order parameter q. The continuous line is a guide for the eyes. 



In the absence of calculations of / and Z when K = 3.3 in the literature (for a explicit 
comparison), in Fig. H]we just present the probability distribution histogram for the order 
parameter q, considering D = 8.605, T = 1.50 and a small lattice of L = 25. From the 
two peaks with very similar heights, we see that the simulations are able to cross the large 
free-energy barriers at the phase coexistence. On the other hand, from an usual Metropolis 
simulations, the system would be trapped in metastable states, evolving to the stable phases 
only after very large MC steps. We should mention that qualitatively our results agree quite 
well with those in Fig. 1 (b) of Ref.- for very similar parameters (actually, there the authors 
use a bigger lattice of L = 32 at T = 1.50, with D = 8.6035). We do not make a direct 
quantitative comparison, e.g., by digitalizing the results in Ref.-, simple because in the 
mentioned figure the scales are not explicitly given. By setting L = 32 (for which the 
eigenvalue method is even better because the greater L), we also have found a balanced 
bimodal distribution at the coexistence, with D 8.6035 as in Ref.-. 

As a final observation, we recall that K = 3.0 has been extensively studied under different 
approaches'^. In particular, it has been shown^ that for this parameter value, the ST is 
an efficient method to study the first-order phase transition for 2D square lattices as small 
as L = 20. 
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IV. THE BELL-LAVIS MODEL 



As a last example, we consider the Bell-Lavis water model^ 3 -* 3 ^. It is defined on a triangular 
lattice, whose occupational variable (i — 1, 2, . . . , V") takes the values 1 (0) if the site is (is 
not) occupied by a water molecule. Moreover, each molecule is described by an orientational 
state, indicating to which nearest neighbor a hydrogen bonding can be formed. So, for a 
non-empty site i, we define the quantity r\\ where j runs over the first neighbors j of i. 
If there is a bonding arm pointed from % towards j, then t\ 3 = 1, otherwise t\ 3 = 0. Two 
adjacent molecules always interact via Van der Waals forces, whereas they do form hydrogen 
bonds provided rf x rj 4 = 1. 

The model is described, in the grand-canonical ensemble, by the Hamiltonian 

U = - a i a 3 ( e hb if + tvdw) - nJ^CTi, (1 1) 

<i,j> i 

where p is the chemical potential and e v d w and are, respectively, the Van der Waals 
and hydrogen bonds interaction energies. The Van der Waals force tends to increase the 
system density by filling the lattice with molecules. On the other hand, the hydrogen bond 
interaction essentially favors an increasing in the hydrogen bonds, so it effectively may limit 
the molecule density if p is negative and small. Finally, the T(Sk) elements are given by 

L 

T(S k ) = exp (a ijk (a ijk + 2a i+1:k )(e vdw + e hb T^ k r i+ i tk + /x)) . (12) 
i=i 

The Hamiltonian in Eq. ( TTT]) exhibits a very rich phase diagram. For instance, recent 
numerical simulations 35 show that the parameter conditions allowing the existence of two 
stable liquid phases take place when the hydrogen bonds are at least three times higher 
than the Van der Wall interaction, or ( = e vdw /e hb < 1/3. In such case, for low negative 
/x values the system presents only a gas phase. By increasing /x a low-density-liquid-phase 
(LDLP) arises. In the limit of higher chemical potentials, we have a high-density-liquid- 
phase (HDLP). At T = 0, the LDLP (HDLP) has a global density of p = 2/3 (p = 1), 
with the hydrogen bond density per molecule being = 3/2 (p^ = 1). Furthermore 3 ^, 
both phase transitions are of first-order: that between the gas and the LDLP occurring 
at p* = -3 (1 + 0/2 and that between the LDLP and HDLP at p* = For T > 0, 

the former remains first-order, ending at a tricritical point, whereas the latter becomes 
second-order, belonging to the Ising universality class^ 5 -. 
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FIG. 5. For the Bell-Lavis model with L = 18, / x p obtained from the integration of the Gibbs- 
Duhem relation (continuous lines) and from the largest eigenvalue of T (symbols). The different 
curves have been offset for a better visualization. 

Next, we will focus on the first-order phase transition case, gas-LDLP, assuming ( = 1/10 
and presenting all the results in units of e^. First, we test the efficiency of the transfer 
matrix to yield the free-energies, hence the weights for the ST. In Fig. |5]we compare some 
values of / calculated from the T approach with those obtained from numerical integration 
of the Gibbs-Duhem equation, S dT — V dp + N dp = 0, at fixed temperatures. Thus, 
dp = pdp where the pressure is related to the free-energy per volume (grand-canonical) by 
the expression / = —p. As it can be seen, even for a small lattice size of L = 18 and low 
T's, the agreement is very good. 

Now, we assume a lattice of L = 24 and a rather small T = 0.25, difficult to simulate 
by standard one-flip algorithms. In Fig. [6] we plot the density p (the order parameter) as 
function of the chemical potential p (the control parameter) around the phase transition 
point. We consider both the ST as well as an usual Metropolis algorithm. Contrary to 
the latter, the ST predicts the phase transition without displaying any hysteresis. In fact, 
hysteresis is a characteristic behavior of methods not able to properly sample the system 
when the phase space presents high free-energy barriers. Also, for the low T considered, 
the transition p* should not be too different from -1.65, the thermodynamic value at T = 
when ( = 1/10. Indeed, an expectation confirmed by Fig. 6. 

Due to the lack of studies discussing first-order phase transition for the Bell-Lavis model 
(either by means of general or dedicated methods), here we compare the proposed ST algo- 
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FIG. 6. For the Bell-Lavis model with L = 24 and at T = 0.25, q versus p calculated from the 
ST (square) and from a common Metropolis method. The hysteresis is due to the latter algorithm 
difficulty in properly sampling the system. 
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FIG. 7. For the Bell-Lavis model with T = 0.25, L = 24 and /i = -1.6533, the probability 
distribution histogram of the order parameter p calculated from the ST and PT methods. The 
continuous line is a guide for the eyes. 

rithm with calculations based on the parallel tempering (PT) approach, recently shown to 
be a very efficient tool to analyze such thermodynamical regime^'^'^'^. In Fig. [7] we plot 
the histogram of the order parameter p at the phase coexistence for T = 0.25 and L = 24. 
In 'tuning' the chemical potential so to have the two peaks of about the same high, we find 
p = —1.6533. We emphasize the quite good agreement between the methods, both being 
able to circumvent metastable states and to promote frequent visits between the gas phase 
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FIG. 8. For the Bell-Lavis model, the order parameter p and compressibility \T versus p for 
different system sizes at T = 0.25. The insets show the collapsed data, respectively, plotted as 
p x (p — p*)V and xt/V x (p _ M*) V. 

and the LDLP. 

We finally perform a finite-size analysis to determine the phase coexistence. We plot in 
Fig. |H]the order parameter p and the compressibility xt as functions of p, assuming different 
system sizes L and fixed T = 0.25. Note that we obtain very smooth curves and that xt 
exhibits sharper peaks as L increases. Moreover, all the isotherms for the density p cross 
each other at the same point p* = —1.6528(1) (as it should be2£, a value close but smaller 
than that for a finite system, e.g., the one in Fig. [7]). After locating p* with accuracy, we 
can perform a rescaling of the data, finding a very good collapse by plotting p x (p — p*) V 
and Xt/V x {f 1 ~ A 4 *) V- Such results confirm the amenability of the present procedure for 
estimating the transition points in discontinuous phase transitions. 

V. REMARKS AND CONCLUSION 

In this paper we have considered an alternative simple protocol to calculate very accurate 
free-energy weights for the ST (simulated tempering) method. It consists in estimating / 
by means of the largest eigenvalue of the transfer matrix—, a quantity directly obtained 
from standard Monte Carlo simulations. To illustrate the approach, we have addressed 
strong first-order phase transitions for different lattice models, namely, Blume-Capel, BEG 
and Bell-Lavis. In such regime, distinct regions of the phase space can be separated by high 
entropic barriers, thus demanding efficient sampling procedures. In all cases, our results have 
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agreed quite well with available precise calculations in the literature, with the advantage that 
our ST has a simple implementation. Next, two remarks are in order. 

First, as already mentioned, an advantage of the present implementation is the relative 
low computation cost in obtaining the refined weights - based on standard MC simulations. 
It is not our purpose here a detailed benchmark comparison between different methods. 
Nevertheless, without going into the merit of the (good) performance of other approaches 
such as Wang-Landau and parallel tempering (for the later, see, e.g., Ref.— ), we observe 
the following. Algorithmically speaking, in principle the ST may be faster since it does not 
demand: either (i) to know the density of states as in the Wang-Landau (usually a hard 
task, e.g., in our examples involving 3 V configurations); or (ii) to simultaneously simulate 
many replicas as in the PT (in the ST only one system realization is considered). 

Second, by increasing too much the system size, the accepting probabilities for the tem- 
perature exchanges decrease, eventually leading to a poor estimation for the thermodynamic 
quantities (but for a possible way to assuage it, see the end of Section II). In fact, the dif- 
ficulty in simulating large systems is not a peculiarity of the ST. It also may be the case 
in others methods like the PT and Wang-Landau. A finite-size analysis circumvent this 
problem, but in practice should use L's up to a certain maximum value L m , for which the 
considered method still works well. The crucial point is whether L m allows a correct ex- 
trapolation to the thermodynamic limit. In our many examples, for L around 25 we already 
can predict this limit. But other situations would require, say, L = 100, which in principle 
could be calculated with our approach if we use a larger set of temperatures {T n }, but with 
AT = (T n — Ti) fixed. Presently, such issue is under investigation and will the subject of a 
forthcoming publication. 

In summary, the examples given show that the ST algorithm allied to our straightforward 
way to calculate the weights is able to deal with free-energy barriers, common at the phase 
coexistence, and therefore can be very useful in characterizing first-order phase transitions. 
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Appendix A: The transfer matrix method to calculate Z 



Here we present a general overview on the transfer matrix method^, and how it can be 
used to calculate Z. 

Consider a general Hamiltonian of a regular lattice, written as 

n = J2nSk,s k+1 ), (Ai) 

k=l 

for Sk = 02,&j • • • > °"L,fc) the state configuration of the k-th layer (each having L sites). 
Also, assume periodic boundary conditions, so that Sk+i = Si- 

From the definition of the (grand-canonical) partition function Z, we have that the prob- 
ability P(S\, S2, Sk) for the layer 1 to have the configuration Si, the layer 2 the configu- 
ration S2, and so on, is given by (/3 = l/(ksT)) 

pfQ Q v exp[-(3U(Si, S 2 )]x...x exp[-(m(S K -i, S K )} x exp[-PH(S K , S ± )} 

• • • , OK) — 7? • 

(A2) 

We then can define the transfer matrix T such that its elements T(S', S") are equal to 
exp[— pH(S', S")], for S' and S" being the configurations of two successive neighbor layers. 
Hence, the r.h.s. of Eq. (|A2]) reads T(S U S 2 ) . . . T(Sjc, Si)/^- 
Next, since 

£ P(S 1 ,...,S K ) = l, (A3) 

Si,...,Sk 

it naturally follows that 

Z = E TOSi, 5 2 ) T(S 2 , S 3 ) ... T{S K -i, S k ) T(S k , Si). (A4) 

Si,...,Sk 

But observe that J2s k 7~(Sk-i, Sk) T(Sk, Sfc+i) = T^(Sk-i, Sfe+i)> with this last term denoting 
the element (Sk-i, Sfc+i) of the matrix T 2 . Thus, using such relation recursively in Eq. (1A4I) . 
one finds 

Z = Y,T k (S 1 ,S 1 ) = Tt[T k }. (A5) 

Si 

If now we calculate all the eigenvalues {A^} of 7~, from basic linear algebra we get 

Z = E(A (m) )^. (A6) 



Finally, for K large enough, the most important contribution in Eq. (1A6I) comes from the 
largest eigenvalue of T, A 1 - ), and thus we recover Eq. ([2]) of Sec. II (which is exact in the 
thermodynamic limit of K — >■ 00). 
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To derive an expression for \(°\ observe that the marginal probabilities 
P(S 1 )= Yl P(S 1 ,S 2 ,S 3 ,...,S N ), P(Sx,S 2 )= ]T P(S 1 ,S 2 ,S 3 ,...,S N ), (A7) 

S2,---,Sn ^3,...,5jv 



are readily obtained from Eq. (1A2|) and from proper products of the T matrix elements, or 
(S 1 = S', S 2 = S") 

P(S') = ' Kb ' b \ P(S', S") = l^> b )l (A8) 

Zj Zj 

For |m) being the eigenvector associated to the eigenvalue A*" 1 - 1 of T, we have that the usual 
spectral expansion of an operator, T = Y, m \m)(m\, yields 

T(S', S") = ]T A (m) ^ m \S') <f) im) *(S"), (A9) 

m 

for cj)( m \S') an element of \m) in the representation of the layer configurations {S}. Therefore 

p(s') = \ }2(* im) ) K <i> {m) (s')<i> imr (s'), ( A1 o) 

Z m 

and 

P(S',S") = yT(.S',5 /, )E( A(m) )^V (m) (^ , )0 (mr (^")- (AH) 

m 

Again, considering K large enough, Eqs. (1A10I) and (1A11I) can be approximated by 

P(S') = 4> {0) (S')^ 0) *(S'), P(S',S") = ^T(S',S")^ \S')4>^\S"). (A12) 
Setting S' = S" in Eq. (1A12[) . we arrive at 

P(S', S') = £ S s ,, s »P(S', S") = t^T(S', S')P(S'). (A13) 



Lastly, summing Eq. ( 1A13I) over S' and identifying the averages (for T(S') = T(S', S')) 



(5s>,s») = E 8w P(S', S"), (T(S')) = £ T(S') P(S'), (A14) 

S'S" S' 



we obtain Eq. ([5]). 
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